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ABSTRACT 

We have modelled the nonlinear development of the secular bar-mode insta- 



bility that is driven by gravitational radiation-reaction (GRR) forces in rotating 
neutron stars. In the absence of any competing viscous effects, an initially uni- 
formly rotating, axisymmetric n = 1/2 polytropic star with a ratio of rotational 
to gravitational potential energy T/|W| = 0.181 is driven by GRR forces to a 
bar-like structure, as predicted by linear theory. The pattern frequency of the 



bar slows to nearly zero, that is, the bar becomes almost stationary as viewed 
from an inertial frame of reference as GRR removes energy and angular momen- 
tum from the star. In this "Dedekind-like" state, rotational energy is stored as 
motion of the fluid in highly noncircular orbits inside the bar. However, in less 
than 10 dynamical times after its formation, the bar loses its initially coherent 
structure as the ordered flow inside the bar is disrupted by what appears to be a 
purely hydrodynamical, short-wavelength, "shearing" type instability. The grav- 
itational waveforms generated by such an event are determined, and an estimate 
of the detectability of these waves is presented. 

Subject headings: neutron stars - - hydrodynamics - - secular instabilities - 
gravitational radiation 



1. Introduction 

As Chandrasekhar (1969) and Tassoul (1978) have discussed in depth, if a star rotates 
sufficiently fast - to a point where the ratio of rotational to gravitational potential energy 
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in the star T/|W| > 0.27 - it will encounter a dynamical instability that will result in the 
deformation of the star into a rapidly spinning, bar-like structure. Although originally iden- 
tified in configurations that are uniformly rotating and uniform in density, more generalized 
analyses made it clear that the dynamical bar-mode instability should arise at approximately 
the same critical value of T/|W| in centrally condensed and differentially rotating stars (cf., 
Ostriker & Bodenheimer 1973). A number of groups have used numerical hydrodynamic 
techniques to follow the nonlinear development of this bar-like structure in the context of 
the evolution of protostellar gas clouds (Tohline, Durisen, & McCollough 1985; Durisen et 
al. 1986; Williams & Tohline 1988; Pickett et al. 1998; Cazes & Tohline 2000) and in the 
context of rapidly rotating neutron stars (New, Centrella & Tohline 2000; Brown 2000). Very 
recently numerical simulations by Centrella et al. (2001) and Shibata, Karino, & Eriguchi 
(2002, 2003) have shown that low-order nonaxisymmetric instabilities can become dynami- 
cally unstable at much lower values of T/|W| in stars that have rather extreme distributions 
of angular momentum. Through linear stability analyses, Karino & Eriguchi (2003) and 
Watts, Andersson & Jones (2004) have attempted to show the connection between these 
instabilities and the classical bar-mode instability discovered in stars with less severe distri- 
butions of angular momentum. In our present analysis we will not be directly investigating 
the onset or development of these dynamical instabilities. 

Classical stability studies have also indicated that a uniformly rotating (or moderately 
differentially rotating) star with T/|W| > 0.14 should encounter a secular instability that 
will tend to deform its structure into a bar-like shape if the star is subjected to a dissipative 
process capable of redistributing angular momentum within its structure. The nonlinear de- 
velopment of this secular instability has not previously been modelled in a fully self-consistent 
manner, so it is not yet clear whether stars that encounter this type of instability will evolve 
to a structure that has a significant bar-like distortion. In this paper, we present results from 
a numerical simulation that has been designed to follow the nonlinear development of the 
secular bar-mode instability in a rapidly rotating neutron star. A force due to gravitational 
radiation-reaction (GRR) serves as the dissipative mechanism that drives the secular devel- 
opment of the bar-mode. By following the development of the bar to a nonlinear amplitude 
and calculating the rate at which angular momentum and energy are lost from the system 
due to gravitational radiation, we are able to provide a quantitative estimate of the dis- 
tance to which such a gravitational-wave source could be detected by existing and planned 
experiments, such as the laser interferometer gravitational- wave observatory (LIGO). 

Chandrasekhar (1970) was the first to discover that gravitational radiation-reaction 
forces can excite the secular bar-mode instability in uniformly rotating, uniform- density 
stars with incompressible equations of state (i.e., the Maclaurin spheroids). This work was 
generalized by Friedman & Schutz (1978) and Comins (1979a,b), to show that the GRR 
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instability extends to stars with any equation of state, and to other nonaxisymmetric modes 
with higher azimuthal mode numbers (m > 2). Managan (1985), Imamura, Friedman, & 
Durisen (1985), Ipser & Lindblom (1990), and Lai & Shapiro (1995) have all shown that the 
critical value of T/\W | at which the GRR secular instability in the (m = 2) bar-mode sets in 
does not depend sensitively on the polytropic index of the equation of state or the differential 
rotation law of the star. These stability analyses have also been generalized to systems in 
which the star is governed by relativistic, rather than purely Newtonian, hydrodynamics 
and gravitational fields (Friedman 1978; Lindblom & Hiscock 1983; Cutler 1991; Cutler & 
Lindblom 1992; Stergioulas & Friedman 1998; Shapiro & Zane 1998; Di Girolamo & Vietri 
2002). 

Lindblom & Detweiler (1977) first showed that viscous processes within compact stars 
can act to suppress the GRR-driven secular bar-mode instability. Various physical viscosities 
have been considered, for example, shear viscosity due to electron and neutron scattering 
(Flowers & Itoh 1976; Cutler & Lindblom 1987), shear viscosity due to neutrino scatter- 
ing (Kazanas & Schramm 1977; Lindblom & Detweiler 1979; Thompson & Duncan 1993), 
bulk viscosity due to weak nuclear interactions (Jones 1971; Sawyer 1989; Ipser & Lindblom 
1991; Yoshida & Eriguchi 1995), or "mutual friction" effects in a superfluid (Lindblom & 
Mendell 1995). In the present work, we will not be investigating the influence of viscous pro- 
cesses on the GRR-driven bar-mode instability, focusing instead on following the nonlinear 
development of the bar-mode in purely inviscid systems. 

Detweiler & Lindblom (1977) and Lai & Shapiro (1995) have made efforts to follow 
the non-linear evolution of rotating neutron stars that are susceptible to the secular (or the 
dynamical) bar-mode instability by using energy and angular momentum conservation to 
construct a sequence of quasi-equilibrium, ellipsoidal configurations. Here, we follow the 
GRR-driven evolution of the bar-mode in an even more realistic way by integrating forward 
in time the coupled set of nonlinear partial differential equations that govern dynamical 
motions in nonrelativistic fluids, and by including a post-Newtonian radiation-reaction force 
term in the equation of motion. Reviews of this and other instabilities that are expected 
to arise in young neutron stars have been written by Lindblom (1997), Lindblom (2001), 
Stergioulas (2003), and Andersson (2003). 



2. Methods 

Using the Hachisu (1986) self-consistent-field technique, we constructed two initial equi- 
librium stellar models governed by Newtonian gravity and an n = 1/2 polytropic equation 
of state; that is, the gas pressure p and density p were related through the expression 
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p = Kp 1+1 ^ n , where K is a constant. Model "SPH" was initially nonrotating and, hence, 
spherically symmetric; model "ROT181" was initially axisymmetric, and uniformly rotating 
with a ratio of rotational to gravitational potential energy T/| W| = 0.181. Other properties 
of these two initial models are detailed in Table 1: M is the mass of the star; R cq and -R po i c 
are the star's equatorial and polar radii, respectively; p is the star's mean density; Vt TOt is the 
angular velocity of rotation; and J is the star's total angular momentum. Columns 2 and 4 
of Table 1 give the values of these various quantities in dimensionless code units, where we 
have assumed the gravitational constant, the star's central density, and the radial extent of 
the computational grid are all equal to unity (i.e., G = p c = G7 gr id = 1). Columns 3 and 5 of 
Table 1 give the value of each quantity in cgs units, assuming both stars have M = 1.4 M Q 
and K = 1.83 x 10~ 10 cm 8 g~ 2 s~ 2 . (This value of K produces a spherical 1.4 M star with 
a radius of 12.5 km, which is characteristic of a neutron star.) 

Each model was introduced into our hydrodynamical code along with a low-amplitude, 
nonaxisymmetric perturbation that was designed to closely approximate the eigenfunction of 
the £ = m = 2 "bar-mode" in a spherical, n = 1/2 polytrope (Ipser & Lindblom 1990). As an 
illustration, Fig. 1 shows the perturbed velocity field that was introduced in the equatorial 
plane of model SPH along with a low-amplitude, bar-like distortion in the density which 
oriented the bar along the vertical axis. Then the nonlinear hydrodynamical evolution of 
each model was followed using the numerical simulation techniques described in detail by 
Motl, Tohline, & Frank (2002). More specifically, we integrated forward in time a finite- 
difference approximation of the following coupled set of partial differential equations: 

^ + V-(H = 0, (1) 

/ dv - , - - , 

P(^T + ^-Vt/) = -Vp-pV($ + K$ GR ), (2) 

^ + v-M = 0, (3) 

V 2 $ = AirGp, (4) 

where v is the fluid velocity, $ is the Newtonian gravitational potential, r = (ep) l / T is 
the entropy tracer, e is the specific internal energy, p =(T — l)ep and T = 1 + 1/n = 3. 
Because the models were initially constructed using a polytropic index n — 1/2 and they 
were evolved using an adiabatic form of the first law of thermodynamics (Eq. 3) with an 
adiabatic exponent T = 3, the models effectively maintained uniform specific entropy at a 
value specified by the initial model's polytropic constant, K. 

In the equation of motion, Eq. (2), we included the post-Newtonian approximation 
to the gravitational radiation-reaction potential produced by a time-varying, £ = m = 2 
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Fig. 1. — Velocity vectors in the equatorial plane of model SPH at time t — 0, but after the 
nonrotating model was perturbed by the £ = m = 2 "bar-mode" eigenfunction drawn from 
linear theory. 
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mass-quadrupole moment (Ipser & Lindblom 1991), 



•» s -V^(i)^^- (5) 

(5) 

where D 22 is the fifth time-derivative of the quadrupole moment and c is the speed of light. 
For modeling purposes, a dimensionless coefficient k was affixed to the radiation-reaction 
potential term in the equation of motion. By adjusting the value of K, we could readily 
remove or artificially enhance the effect of this non-Newtonian GRR force. 

As implemented on our cylindrical computational mesh (w, <p, z), D 22 and its first time- 
derivative were evaluated using the following expressions, 

D 22 = \[^jp^~^d 3 x, (6) 

D% = \hr i Mvw-iv^\e-^d\. (7) 

V »7T J 

Following Lindblom, Tohline, & Vallisneri (2001, 2002) we have assumed that the quadrupole 
moment has a time-dependence of the form D 22 oc e~^ 22 ', hence, 

D$ = (-ioo 22 ) n D 22 , (8) 

where, at any instant in time, the complex eigenmode frequency uj 22 = uj v + iuj\ can be 
determined by taking the ratio D 22 /D 22 , so in Eq. (5) we set 

D$ = \uto\*D$. (9) 



3. Predictions of Linear Theory 

Using the linear perturbation techniques described by Ipser & Lindblom (1990, 1991) 
we determined that in model SPH, u r = Re(cc>22) = ±(1.21 ± 0.01)f2 an d, if the model is 
scaled to a mass of 1.4M and a radius of 12.5 km, ui = lm.(u 22 ) = —1.00 ± 0.01 x 10~ 3 k,Q , 
where Qq = \pnGp = 1.308 in dimensionless code units. (The uncertainty is estimated 
from the values that are determined numerically by the linear perturbation method for 
different radial resolutions.) In model SPH, therefore, we should expect the amplitude of 
the bar-mode to damp exponentially on a time scale r GR = 1/1^1 = 193K _1 r pattcrn , where 
^pattern = 27r/|c<j r | = 3.97 in dimensionless code units. 

Linear perturbation analyses have not yet provided quantitative values of the bar-mode 
eigenfrequency in rapidly rotating, n — 1/2 polytropes. From the information given in Ipser 
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Fig. 2. — Solid curves depict the time-evolution of the amplitude I-D22I (bottom) and the real 
(top) and imaginary (middle) components of the £ = m = 2 bar-mode frequency from model 
SPH. Time is shown in units of the predicted pattern period; I-D22I has been normalized 
to MRq • and frequencies are shown in dimensionless code units. Dash-dotted lines show 
predictions of linear theory. 



& Lindblom (1990, 1991), however, we expect that in model ROT181, (a) uj t /Qq should be 
positive but close to zero, as viewed from an inertial reference frame; and (b) Ui/Qq should 
be slightly positive - that is, the mass-quadrupole moment should grow exponentially, but on 
a time scale that is very long compared to the damping time predicted for the nonrotating 
model, SPH. More specifically, if |w r |/r2o proves to be an order of magnitude smaller in 
model ROT181 than it is in model SPH, then we should expect tqr to be ~ 10 5 times larger 
because the amplitude of the GRR driving term in Eq. (2) is proportional to ui\ 2 - 



4. Barmode Evolutions 

4.1. Model SPH 

Initially the perturbation applied to model SPH had an amplitude Dq = \D 2 2{t = 0)| ~ 
10~ 3 and a velocity field (see Fig. 1) designed to excite the "backward moving" £ = m = 2 
bar-mode, that is, the mode for which uu r < 0. We followed the evolution of the model on a 
cylindrical grid with a resolution of 66 x 128 x 130 zones in w, <p, and z, respectively, and 
with the coefficient of the radiation-reaction force term in Eq. (2) set to the value k = 20. 
The effect of this was to shorten the timescale for the exponential decay by a factor of 20, to 
a predicted value tqr = 9.63r pat tcrn- By shortening the decay timescale in this manner, we 
were able to significantly reduce the amount of computational resources that were required 
to follow the decay of the barmode while maintaining a decay rate that was slow compared to 
the characteristic dynamical time of the system, Qq 1 ?» 0.2r pat tcrn- This is the same technique 
that was successfully employed by Lindblom et al. (2001, 2002) in an earlier investigation of 
the r-mode instability in young neutron stars. 

Figure 2 displays the key results from our SPH model evolution. The solid curves in the 
top two frames display u r and u)\ as a function of time through just over nine pattern periods. 
Each of these frequencies oscillate about a fairly well-defined, mean value: (uj r ) ss —1.56 = 
— 1.190o; (ui) ~ —0.03 = — 0.023f2o. Oscillations about these mean values initially had an 
amplitude ~ ±0.05 = 0.038f2o, indicating that our initial nonaxisymmetric perturbation did 
not excite a pure eigenmode, but these oscillations decreased in amplitude somewhat as the 
evolution proceeded. Our measured values of uj r and u>{ are within 3% and 15%, respectively, 
of the values predicted from linear theory (see §3). The solid curve in the bottom frame of 
Fig. 2 shows in a semi-log plot the behavior of I-D22I with time. (Note that in this figure, 
I-D22I has been normalized to MR^ q .) There is a clear exponential decay with a measured 
damping time (given by the slope of the solid curve) t G r ~ 8.45r pattern . This decay time 
is completely consistent with the measured value of (u[) that we have obtained from the 
middle frame of Fig. 2 and, again, within 15% of the predicted value (illustrated by the solid 
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dash-dotted line in the top frame of the figure). The somewhat larger discrepancy in the 
measured value of uji is most probably due to the fact that the GRR formalism used here 
was derived under the assumption that \u>i\ << \uj r \. Since u>i is caused by the $gr potential 
which is proportional to the fifth power of the frequency, fractional discrepancies which are 
of order 5\uii/ui r \ ps 0.1 are not unexpected. 

The first row of numbers in Table 2 summarizes these simulation results. Specifically, 
columns 4, 5 and 6 list the values of u r , u>i, and tqr that have been drawn directly from 
Fig. 2; all three numbers are given in dimensionless code units. In the last two columns of 
this table, the real and imaginary frequencies have been reexpressed in units of the dynamical 
frequency, Qq. In the last column, we also have adjusted Ui by the factor of k in order to 
show the frequency (and associated growth rate) as it would appear in a real neutron star 
where the GRR force would not be artificially exaggerated. 



4.2. Model ROT181 

4-2.1. Radiation-reaction with K = 1.75 x 10 5 

Model ROT181 was introduced into our hydrocode with a nonaxisymmetric perturbation 
in the density that had the same structure as the perturbation that was introduced into model 
SPH. Because we expected the natural oscillation frequency of the bar-mode to be close to 
zero (as viewed from an inertial reference frame), however, we did not perturb the velocity 
field of the model. We followed the evolution of model ROT181 on a cylindrical grid with a 
resolution of 130 x 128 x 98 zones in zu, <f), and z, respectively, and with the coefficient of the 
radiation-reaction force term set to k — 1.75 x 10 5 . Note that fewer vertical grid zones were 
required than in model SPH because model ROT181 was significantly rotationally flattened, 
but more radial zones were used than in model SPH in order to allow room for model 
ROT181 to expand radially during the nonlinear-amplitude phase of its evolution. A much 
larger value of k was selected because, as explained earlier, the natural growth rate of the 
bar- mode in model ROT181 was expected to be orders of magnitude smaller than the decay 
rate measured in model SPH. 

Figures 3 and 4 display some of the key results from this ROT181 model evolution. The 
bottom frame of Fig. 3 shows the time-dependent behavior of the real (dash-dotted curve) 
and imaginary (solid curve) components of u>22, in our code's dimensionless frequency units; 
the solid curve in the top frame displays the time- dependent behavior of I-D22I, normalized 
to MR^ . Figure 4 shows how the global parameters T/|VF| (solid curve) and J (dashed 
curve) evolve with time. The behavior of the model can be best described in the context of 



-10- 



Table 1. Initial Model Parameters 





Model SPH 


Model ROT181 


Units 


code 


cgs 


code 


cgs 


c 


3.122 


3.00 x 10 10 


1.634 


3.00 x 10 10 


M 


1.359 


2.80 x 10 33 


0.1716 


2.80 x 10 33 


K 


0.7825 


1.83 x 10~ 10 


0.1281 


1.83 x 10~ 10 


-rteq 


0.8413 


1.25 x 10 6 


0.6102 


1.97 x 10 6 


-'tpolc 


0.8413 


1.25 x 10 6 


0.2746 


8.86 x 10 5 


P 


0.5460 


3.42 x 10 14 


0.4922 


2.39 x 10 14 


J L ro i 


0.0 


0.0 


0.9705 


5.52 x 10 3 


J 


0.0 


0.0 


0.01632 


1.58 x 10 49 



Table 2. Simulation Results 



Model 



Q 



L) r 



0Ji r GR uj r /tt u-J^qk) 



SPH 1.308 

ROT181 1.488 



2.00 x 10 1 
1.75 x 10 5 



-1.56 
+0.27 



-0.03 
+0.08 



34 
12 



-1.19 
+0.18 



-1.1 x 10~ 3 
-3.1 x lO" 7 
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Fig. 3. — The time-evolution of the amplitude I-D22I (top) and the £ = m = 2 bar- mode 
frequency CJ22 (bottom) from two rapidly rotating models. Time is shown in units of the 
initial rotation period T sp i n = 27r/O rot of the model; I-D22I has been normalized to Mi?^ q ; and 
frequencies are shown in dimensionless code units. Curves that terminate at approximately 
17r sp i n display data from model ROT181 and curves that extend past 35r spin show data from 
the lower resolution model ROT 179. In the bottom panel, both the real (eg., dash-dotted 
curve for model ROT181) and imaginary (eg., solid curve for model ROT181) components 
of co>22 are displayed. 
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three different evolutionary phases: Early [0 < t/r sp i n < 7]; intermediate [7 < t/r sp i n < 12]; 
and /ate [t/r sp i n > 12], where r spin = 2ir/Q rot = 6.47 in dimensionless code units. 

During the model's early evolution, both components of the frequency ^22 oscillate about 
well-defined, mean values: (u T ) ss 0.27 = 0.181O ; (u>i) ~ 0.08 = 0.054fi . (Following Ipser 
& Lindblom 1991, we define Qq in terms of the mean density po of a spherical star that has 
the same M and K as model ROT181, that is, VLq = \fnGpQ = 1.488 in dimensionless code 
units; see Table 2.) During this same phase of the evolution, both J and T/|W| remain fairly 
constant, but \D 2 2\ increases exponentially with a growth time (obtained from the slope of 
the displayed curve) tgr ~ 1.85r spin . This growth time is completely consistent with the 
measured value of (ui), from which we would expect TGR/r spin = (ui)~ 1 (Q rot /2n) = 1.93. 
The second row of numbers in Table 2 summarizes these simulation results. 

After approximately seven rotation periods, the amplitude of I-D22I begins to saturate, 
and the model deforms into a clearly visible bar-like configuration with an axis ratio measured 
in the equatorial plane of approximately 2:1 (see Fig. 5). The bar-like structure is initially 
spinning with a frequency given by (u T )/2, as measured during the early phase of the ROT181 
evolution. This pattern frequency of the bar is a factor of 7.2 smaller than the rotation 
frequency Q vot of the model in its initial, axisymmetric state, so it is not surprising that 
the bar also exhibits sizeable internal motions - it has a "Dedekind-like" structure. Figure 
5 illustrates the structure of the model at this time. Both frames contain the same set of 
equatorial-plane, isodensity contours delineating the bar, along with a set of velocity vectors 
depicting the fluid flow inside the bar: on the left-hand-side, the velocity vectors are drawn 
in a frame corotating with the bar (i.e., rotating at the frequency, (u T )/2) to illustrate the 
elliptical streamlines of fluid flow within the "Dedekind-like" bar; on the right- hand- side, 
the velocity vectors are drawn in a frame rotating at the frequency Jl rot . When viewed in 
this latter frame, one sees a global velocity structure that is very similar to the flow-field 
depicted in Fig. 1, that is, it resembles the natural eigenfunction of the £ = m = 2 bar- mode 
that was derived by perturbation analysis for nonrotating spherical stars, such as our model 
SPH. We note that this velocity structure developed spontaneously in model ROT181, as 
the initial model contained no velocity perturbation. 

During this intermediate phase of the model's evolution, the bar remains a robust con- 
figuration, but its pattern frequency slows as the system loses approximately 10% of its 
angular momentum (through gravitational radiation) and T/|W| drops to a value ~ 0.156. 
It is particularly interesting to note that, during this phase of the evolution, the GRR driv- 
ing term in the equation of motion reaches a maximum, then drops as rapidly as it initially 
rose; this is illustrated in Fig. 6, where we have plotted the time-dependent behavior of the 
product, |ct-?22 1 5 1 -D22 1 - Although the bar maintains a nonlinear structure, i.e., I-D22I remains 
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large, during this intermediate phase of the model's evolution $gr drops quickly in concert 
with a decrease in the frequency (u^l- 

During the late phase of the model ROT181 evolution, the Dedekind-like bar began 
to lose its coherent structure. Small-scale fluctuations in the density and velocity fields 
developed throughout the volume of the bar, and these fluctuations grew in amplitude on a 
dynamical time scale. Even vertical oscillations developed throughout the model, disrupting 
both the vertically stratified planar flow and reflection symmetry through the equatorial 
plane that persisted throughout the early and intermediate phases of the model's evolution. 
After approximately 15r sp i n , the model was no longer a recognizable bar, although it remained 
decidedly nonaxisymmetric, showing density and velocity structure on a wide range of scales 
in all three dimensions. Figure 7 provides a snapshot of model ROT181's structure at 
t = 19.9 T spin during the late phase of its evolution. (Actually, Fig. 7 is drawn from the late 
phase of a "revised" evolution of model ROT181, which was evolved further in time; see 
§4.2.3 for details.) Isodensity contours reveal a nonaxisymmetric structure that no longer 
can be described simply as a bar and, when viewed from a frame rotating at a frequency 
Q rot (the right-hand frame), the flow field is seen to be more complex than in the bar. 



4-2.2. Detectability of gravitational-wave radiation 

A rapidly spinning neutron star located in our Galaxy (and perhaps anywhere in our 
local group of galaxies) that acquires the type of nonlinear-amplitude, bar-like structure 
that developed in model ROT181 will produce gravitational radiation at a frequency and 
amplitude that should soon be detectable by gravitational- wave detectors such as LIGO 
(Abramovici et al. 1992; Abbott et al. 2004), VIRGO (Acernese et al. 2002), GEO600 (Willke 
et al. 2002; Gossler et al. 2002), or TAMA300 (Tagoshi et al. 2001). As our simulation shows, 
however, both the amplitude and pattern frequency of the bar — and, hence, the strength and 
observed frequency of the gravitational radiation — will vary with time. To illustrate this, 
Fig. 8 depicts the evolution of model ROT181 across a "strain- frequency" diagram, which 
is often referenced by the experimental relativity community when discussing detectable 
sources of gravitational radiation. Specifically, the dimensionless strain h norm = \/h 2 + + h\, 
where h + and h x are the two polarization states of gravitational waves. For an observer 
located a distance r along the axis (6 = 0, <fi = 0) of a spherical coordinate system with 
the origin located at the center of mass of the system, we have h + = ^-(l X x ~ lyy) and 
hx = %\$xy, where the reduced moment of intertia l\ m = J p(xix m — ^5i m XkXk)dx 3 . To 
obtain the strain values h norm shown in Fig. 8, we have assumed r = lOkpc, and the time- 
derivative of each reduced moment of inertia was evaluated numerically using the method 
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Fig. 4. — The time-evolution of the angular momentum J and the energy ratio T/|W| 
from model ROT181; J is in dimensionless code units, time is shown in units of the initial 
rotation period of the model. During the intermediate phase of the evolution, both quantities 
noticeably drop as angular momentum is lost via the GRR force term in the equation of 
motion. 
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recommended by Finn & Evans (1990). Model ROT181's evolutionary trajectory in this 
diagram is strikingly similar to the trajectory that was predicted by Lai & Shapiro (1995) 
- see their Fig. 4 - using a much simpler, approximate model for the development of the 
secular bar-mode instability in young neutron stars. 

In order to estimate the distance to which a gravitational wave source of this type 
would be detectable by a gravitational- wave interferometer, such as LIGO, we could integrate 
under the curve in Fig. 8, taking into account the amount of time that the source spends in 
each frequency band. Because we have artificially amplified the strength of the GRR force, 
however, our model evolves through frequency space along the curve shown in Fig. 8 much 
more rapidly than would be expected for a real neutron star that experiences this type of 
instability, hence our model cannot be used directly to estimate the length of time that such a 
source would spend near each frequency. However, Owen & Lindblom (2002) have outlined a 
method by which the detectability of a source can be estimated from a knowledge of A J, the 
total angular momentum that is radiated away from the source via gravitational radiation. 
Specifically, the signal-to-noise ratio S/N that could be achieved by optimal filtering can be 
estimated from the expression, 

SV 4G |AJ| 



N J ^rrmc^r 2 fSh(f) 

where m is the azimuthal quantum number (m = 2 for the bar-mode), r is the distance to the 
source, and Sh(f) is the power spectral density of the detector noise at frequency /. From our 
model ROT181 evolution, we find A J = 1.67 x 10 48 g cm 2 s" ! ; and when it reaches its design 
sensitivity, LIGO's 4 km interferometer noise curve 1 should exhibit \^Sh ~3x 10~ 23 Hz -1 ' 2 
at / = 220 Hz, which is the characteristic frequency of the spinning bar in model ROT181. 
From expression (10), we therefore estimate that a source of the type we are modelling will 
be detectable by LIGO with a S/N > 8 out to a distance of 2 Mpc. (During LIGO's S3 
science run in late 2003, the 4 km LHO interferometer had already come within a factor 
of two of this design sensitivity 1 ) With Advanced LIGO (using sapphire test masses, the 
projected noise curve 2 gives \fS~h ~ 2.0 x 10~ 24 Hz^ 1 ' 2 at / = 220Hz) we estimate that this 
type of source will be detectable with S/N > 8 out to 32 Mpc. 



1 The projected noise curve for LIGO's 4 km interferometers - published as part of the LIGO Sci- 
ence Requirements Document (SRD) - and the actual noise curve achieved by the 4 km interfer- 
ometer at the LIGO-Hanford Observatory (LHO) during the S3 science run can be obtained from 
www.ligo.caltech.edu/~lazz/distribution/LSCJ3ata/. 

2 Projected noise curves for the Advanced LIGO design using either sapphire or silica test masses can be 
obtained from www.ligo.caltech.edu/advLIGO. 
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Of course the detectability of gravitational waves generated by the secular bar-mode 
instability will also depend on the frequency with which such events occur nearby. To 
estimate an event rate we can draw on the discussion of Kokkotas (2004) where an estimate 
was made of the event rate of the dynamical bar-mode instability in young neutron stars. 
Since the conditions required for the onset of the secular bar-mode instability (T/|W| > 0.14) 
are almost as extreme as the conditions required for the onset of the dynamical bar-mode 
instability (T/|W| > 0.27), it would be very surprising if the two event rates were not similar. 
If we assume that only young neutron stars can be rotating rapidly enough to be susceptible 
to either bar-mode instability, and if we assume that a neutron star can form only from the 
collapse of the core of massive star, then a reasonable upper limit on the rate of these events 
will be given by the event rate of Type II supernovae, that is, 1-2 per century per gas-rich 
galaxy (Cappellaro, Evans, & Turatto 1999). (Another scenario is that rapidly rotating 
neutron stars form from the accretion-induced collapse of white dwarfs. But according to 
Liu 2002, the frequency of such events is orders of magnitude lower than the event rate of 
Type II supernovae.) Adopting a local galaxy density of n g ~ 0.01 Mpc~ 3 (Kalogera et al. 
2001), we should expect < 30 Type II supernovae each year out to 32 Mpc. Now not all Type 
II supernovae will produce neutron stars (Kokkotas 2004 estimates, for example that 5-40% 
of supernova events produce black holes instead), and only a fraction / rot of neutron stars 
will be formed with sufficient rotational energy to be susceptible to a bar-mode instability, so 
the predicted event rate should be reduced accordingly. A naive estimation based on angular 
momentum conservation during core collapse suggests that virtually all newly born neutron 
stars will be formed rapidly rotating and, therefore, / rot ~ 1; this is the direction Kokkotas 
(2004) leans. But models of axisymmetric core collapse (Tohline 1984; Dimmelmeier, Font, & 
Miiller 2002a,b; Ott et al. 2004) indicate that the ratio of energies T/| W| in a newly formed 
neutron star is quite sensitive to the equation of state of the core during its collapse and it 
is easy to imagine physical scenarios in which appropriately rapidly rotating neutron stars 
will rarely be formed; therefore, f TOt <C 1. At the present time it is not clear which picture 
is more correct, but adopting the more optimistic view it should be possible for LIGO to 
detect on the order of ten such events each year. 



4-2.3. Model Convergence 

In an effort to determine whether the Dedekind-like bar structure was destroyed during 
the late phase of the ROT181 model evolution as a result of physically realistic, hydrody- 
namical processes, or by a radiation-reaction force that was artificially too large, we set 
K = then re-ran the last segment of the simulation, starting from t = llr sp i n . This "re- 
vised" evolution produced results that were qualitatively identical to the late phase of the 
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GRR-driven evolution. That is, the bar was destroyed by the dynamical development of 
velocity and density structure on a wide range of scales in all three dimensions. In an effort 
to quantitatively describe this relatively complex structure, Fig. 9 shows a representation of 
the azimuthal Fourier-mode amplitudes of the model's density distribution at two points in 
time: t = 10r spin , when the bar was well-developed; and t = 20r spin , after the higher-order 
nonaxisymmetric structure was well-developed. (Note that the late phase of this "revised" 
evolution was followed somewhat farther in time than the original model ROT181 evolution 
described in §4.2.1.) At the earlier time, only the m = 2 amplitude contained a significant 
amount of power, and all odd amplitudes were smaller than their even neighbors. At the 
later time the Fourier-mode amplitudes appear to be related to one another by a simple 
power law, indicating that power has been spread smoothly over all resolvable length scales. 

As an additional test of the reliability of our results, we repeated our rotating model 
evolution on a computational grid that had a factor of two poorer resolution in each of 
the three spatial dimensions; specifically, the new simulation was performed on a grid with 
66x64x66 zones in w, <f), and z, respectively. On this lower resolution grid, it was not possible 
to begin the evolution from precisely the same initial state as model ROT181. However, we 
were able to construct a uniformly rotating, n = 1/2 polytrope with T/|W| = 0.179 (only 
1% less than the corresponding initial energy ratio of model ROT181); we introduced a 
nonaxisymmetric density perturbation that produced approximately the same initial mass- 
quadrupole moment amplitude I-D22I as in model ROT181; and throughout the evolution the 
coefficient of the radiation-reaction force term in Eq. 2 was set to k = 1.75 x 10 5 , as in model 
ROT181. Hereafter, we will refer to this lower resolution evolution as model ROT179. 

The key results of this lower resolution evolution are illustrated by the curves in Figs. 3 
and 6 that extend beyond r spin = 35. As shown in the bottom frame of Fig. 3, the real 
and imaginary components of the eigenfrequency W22 identified during model ROT179's 
early evolution were nearly identical to the values measured in model ROT181; as shown 
in the top of Fig. 3, I-D22I grew exponentially up to a nonlinear value, levelled off, then 
slowly decayed as the bar-mode's pattern frequency slowed; and as shown in Fig. 6 the 
strength of the GRR force grew steadily up to a maximum value that was somewhat lower 
in amplitude and somewhat delayed in time compared to model ROT181, then almost as 
rapidly it dropped by an order of magnitude, as in model ROT181. Finally, during the 
late phase of model ROT179's evolution, the bar's coherent structure was destroyed by the 
development of dynamical structure on much smaller scales, just as was observed during the 
late phase of model ROT181's evolution. (This phenomenon is evidenced in Fig. 3 by the 
rapid oscillations in both components of U22 at times t > 32r sp i n .) In this lower resolution 
evolution, however, the small-scale dynamical structure took roughly twice as long to develop 
as in model ROT181. The somewhat slower initial growth of the bar-mode and the bar's lower 
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peak nonlinear amplitude can both be attributed directly to this simulation's coarser spatial 
resolution. The delay in development of the smaller scale structure was almost certainly due, 
in part, to our inability to resolve structure on the smallest scales in model ROT179, but 
the delay may also have occurred, in part, because the bar, itself, was never as pronounced 
as in model ROT181. Similar behavior has been observed in simulations that have analyzed 
the long-term stability of r- mode oscillations in young neutron stars (Gressman et al. 2002). 



5. Summary and Conclusions 

Using nonrelativistic, numerical hydrodyamical techniques coupled with a post-Newtonian 
treatment of GRR forces, we have simulated the nonlinear development of the secular bar- 
mode instability in a rapidly rotating neutron star. In each simulation we have artificially 
enhanced the strength of the GRR force term in the equation of motion (by selecting values 
of the parameter k > 1) in order to be able to follow the secular development of the bar with 
a reasonable amount of computing resources. In each case however, k was set to a small 
enough value that the amplitude of the mass-quadrupole moment changed slowly, compared 
to the dynamical time scale of the system, thus ensuring that the system as a whole re- 
mained in dynamical equilibrium. We first tested our simulation technique by studying the 
evolution of the I = m = 2 bar- mode in a nonrotating neutron star model (model SPH). The 
developing bar-mode exhibited an azimuthal oscillation frequency within 3% of the frequency 
predicted by linear theory, and the amplitude of the bar-mode damped (as predicted) at a 
rate that was within 15% of the rate predicted by linear theory. 

Next, we evolved a rapidly rotating model (model ROT181), which was predicted by 
linear theory to be unstable toward the growth of the bar-mode. From the early "linear- 
amplitude" phase of this model's evolution, we measured the bar-mode's azimuthal oscil- 
lation frequency and its exponential growth rate; the values are summarized in Table 2. 
The oscillation frequency (uj t )/Q was almost an order of magnitude smaller than in model 
SPH, and (uj[)/(Q K,) was four orders of magnitude smaller than (and had the opposite sign 
of) the value measured in model SPH. Both of these frequency values reflect the fact that 
model ROT181 was rotating only slightly faster than the marginally unstable model (pre- 
dicted to have T/|VF| ~ 0.14), in which both components of c^22 should be precisely zero. 
We watched the unstable bar-mode grow up to and saturate at a sufficiently large, nonlin- 
ear amplitude that the bar-like distortion was clearly visible in two- and three-dimensional 
plots of isodensity surfaces. This nonlinear bar-like structure persisted for several rotation 
periods and, during this intermediate phase of the ROT 181 model evolution, we tracked 
the frequency and amplitude of the gravitational radiation that should be emitted from the 
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configuration due to its time-varying mass-quadrupole moment. Our model's evolution in a 
"strain-frequency" diagram closely matches the evolutionary trajectory predicted by Lai & 
Shapiro (1995), lending additional credibility to their relatively simple (and inexpensive) way 
of predicting the evolution of such systems as well as to our first attempt to model such an 
evolution using nonlinear hydrodynamical techniques. During the late phase of our model 
ROT181 evolution, the bar lost its coherent structure and the system evolved to a much 
more complex nonaxisymmetric configuration. The general features of this late phase of the 
evolution were reproduced when the simulation was re-run on a coarser computational grid, 
and even when the GRR forces were turned off. So while the size and shape of the interme- 
diate phase "Dedekind-like" structure of our model may well have been influenced strongly 
by the excessive strength of the GRR force used in our simulation, it appears as though the 
final complex "turbulent" phase of the evolution was governed by purely hydrodynamical 
phenomena. 

It is not clear what physical mechanism was responsible for the development of the 
small-scale structure and subsequent destruction of the bar during the late phase of the evo- 
lution of model ROT181. Because the bar's structure was "Dedekind-like" - that is, fluid 
inside the bar was moving along elliptical streamlines with a mean frequency that was signif- 
icantly higher than the bar pattern frequency - it is tempting to suggest that the small-scale 
structure arose due to differential shear. But, according to Hawley, Balbus & Winters (1999), 
coriolis forces are able to stabilize differentially rotating, astrophysical flows against shearing 
instabilities even in accretion disks where the shear is much stronger than in our "Dedekind- 
like" bar. (See, however, Longaretti 2002 for an opposing argument.) Furthermore, other 
models of differentially rotating astrophysical bars (Cazes & Tohline 2000; New, Centrella 
& Tohline 2000) do not appear to be susceptible to the dynamical instability that destroyed 
the bar in our ROT181 model evolution. We suspect, instead, that the late-time behavior 
of model ROT181 results either from nonlinear coupling of various oscillatory modes within 
the star, or from an "elliptic flow" instability similar to the one identified in laboratory fluids 
that are forced to flow along elliptical streamlines. The dissipative effect of mode-mode (ac- 
tually, three-mode) coupling has been examined in depth by Schenk et al. (2002) and Arras 
et al. (2003) in the context of the r-mode instability in young neutron stars, and Brink, 
Teukolsky, & Wasserman (2004) have shown the connection between this process and the 
rapid decay of the r-mode in extended numerical evolutions such as the ones performed by 
Gressman et al. (2002). However, this phenomenon has not yet been studied to the same 
degree in relation to the £ = m = 2 f-mode. Lifschitz & Lebovitz (1993), Lebovitz & Lif- 
schitz (1996), and Lebovitz & Saldanha (1999) have demonstrated that the "elliptic flow" 
instability seen in laboratory fluids is likely to arise in self-gravitating ellipsoidal figures of 
equilibrium, especially if they have "Dedekind-like" internal flows. Additional analysis and, 
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very likely, additional nonlinear simulations will be required before we are able to determine 
which (if either) of these mechanisms was responsible for the destruction of the bar in our 
ROT181 model evolution. 

Our nonlinear simulation of model ROT181 demonstrates that when a rapidly rotating 
neutron star becomes unstable to the secular bar-mode instability, the bar-like distortion 
may grow to nonlinear amplitude and thereby become a strong source of gravitational radia- 
tion. However, it will not be a long-lived continuous-wave source, as one might optimistically 
have expected; in our simulation, the nonlinear-amplitude bar survived fewer than ten rota- 
tion periods. In a real neutron star the GRR forces will be much weaker than those of our 
simulation, so we expect the bar mode to grow and persist for many more rotation periods. 
However, we also expect the amplitude of the bar mode to saturate at a much lower am- 
plitude in a real neutron star. Nevertheless, we expect the bar mode to persist in rapidly 
rotating neutron stars long enough to allow gravitational radiation to remove sufficient angu- 
lar momentum for them to relax into a secularly stable equilibrium state. Thus the amount 
of angular momentum radiated away in real neutron stars should be comparable to that in 
our simulation. While such astrophysical systems may not be the easiest sources to detect 
with broadband, gravitational-wave detectors such as LIGO because the frequency of the 
emitted radiation will change steadily with time, our estimates suggest that gravitational 
waves arising from the excited secular bar-mode instability in rapidly rotating neutron stars 
could well be detectable in the not too distant future from neutron stars as far away as 
32 Mpc. 

After submitting this paper for publication, we became aware that Shibata & Karino 
(2004) have just completed an investigation similar to the one presented here in which they 
have utilized post-Newtonian simulations to study the nonlinear development of the secular 
bar-mode instability in rapidly rotating neutron stars. Their initial models were differentially 
rotating, n = 1 (Y = 2) polytropes with 0.2 < T/|W| < 0.26. The early and intermediate 
phases of their model evolutions agree well with the results of our model ROT181 evolution, 
that is, the bar-mode grew exponentially at rates consistent with the predictions of linear 
theory and reached a nonlinear amplitude, producing an ellipsoidal star of moderately large 
ellipticity. The strength of the GRR force used in our simulations was considerably larger 
than theirs. This may explain why the bar mode grows to a larger amplitude and why, in 
turn, there is a more significant decrease in the pattern frequency of the bar as it evolves 
toward a Dedekind-like configuration in our simulation. This may also explain why the bar 
mode structure was ultimately destroyed by short wavelength disturbances in our evolutions 
while such turbulence had not yet developed in theirs. 

We thank Gabriela Gonzalez, Peter Saulson, Peter Fritschel and Kip Thorne for guidance 
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in obtaining the LIGO noise figures used in our analysis. We also thank an anonymous referee 
for recommending several ways in which the presentation of our results could be improved. 
This work was supported in part by NSF grants AST-9987344, AST-0407070, PHY-0326311 
and NASA grant NAG5-13430 at LSU; and by NSF grants PHY-0099568, PHY-0244906 and 
NASA grants NAG5-10707, NAG5-12834 at Caltech. Most of the simulations were carried 
out on SuperMike and SuperHelix at LSU, which are facilities operated by the Center for 
Computation and Technology whose funding largely comes through appropriations by the 
Louisiana state legislature. 
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Fig. 5. — The structure of model ROT181 is shown at time t = 8r spin , during the intermediate 
phase of its evolution. In both frames, solid curves are isodensity contours in the equatorial 
plane while vectors illustrate the equatorial-plane, velocity flow field as viewed from a frame 
rotating with a specific frequency as follows: fif rame = (u T )/2 (left); fi framc = Q rot (right). 
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Fig. 6. — From model ROT181 (ROT179), the solid (dotted) curve depicts the time-evolution 
of the product ^22 1 ^22 1 , which indicates the strength of $gr hi the equation of motion. Time 
is shown in units of the initial rotation period. 
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Fig. 7. — The neutron star's structure is shown at time t = 19.9 r sp i n during the late 
phase of the "revised" ROT181 model evolution. In both frames, solid curves are isodensity 
contours in the equatorial plane while vectors illustrate the equatorial-plane, velocity flow 
field as viewed from a frame rotating with a specific frequency as follows: Of ramc = (left); 

^framc = ^rot (right). 
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Fig. 8. — The solid curve traces the evolution of model ROT181 in a "strain-frequency" 
diagram from 6r spin to llr spin . As is schematically illustrated by the vertical dotted line, 
initially, the amplitude h novm of the gravitational wave signal grows at a constant frequency, 
/ = u) T /(2i{) ps 240 Hz. As energy and angular momentum are radiated from the system, the 
frequency drops monotonically, and the strain reaches a maximum amplitude then steadily 
declines. 
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Fig. 9. — A spectrum of the Fourier-mode amplitude of the azimuthal density distribution 
is shown at time t = 10r spin (filled circles), when the bar was well-developed, and at time 
t = 20r spin (open circles), after the higher-order modes destroyed the coherent bar in the 
"revised" evolution of model ROT 181. To guide the eye, amplitudes determined for various 
modes at the same time are connected by straight line segments. 



